##########################################################################################

library(data.table)
library(optparse)
library(dplyr)
library(RColorBrewer)
library(ggpubr)
library(ggsci)

##########################################################################################

option_list <- list(
    make_option(c("--ccf_file"), type = "character") ,
    make_option(c("--type"), type = "character") ,
    make_option(c("--gene_list"), type = "character") ,
    make_option(c("--sample_info"), type = "character") ,
    make_option(c("--out_path"), type = "character")
)

if(1!=1){
    
    ccf_file <- "~/20220915_gastric_multiple/dna_combinePublic/mutationTime/result/All_CCF_mutTime.addShare.tsv"
    sample_info <- "~/20220915_gastric_multiple/dna_combinePublic/config/tumor_normal.class.list"
    gene_list <- "~/20220915_gastric_multiple/dna_combinePublic/mutsig_check/smg.list"
    out_path <- "~/20220915_gastric_multiple/dna_combinePublic/images/GeneCCF"
    type <- "All"

}

###########################################################################################

parseobj <- OptionParser(option_list=option_list, usage = "usage: Rscript %prog [options]")
opt <- parse_args(parseobj)
print(opt)

gene_list <- opt$gene_list
sample_info <- opt$sample_info
out_path <- opt$out_path
ccf_file <- opt$ccf_file
type <- opt$type

###########################################################################################

dir.create(out_path , recursive = T)
col <- c( "#006699","#DDA520"  )

###########################################################################################

dat_info <- data.frame(fread(sample_info))
dat_gene <- data.frame(fread(gene_list))
dat_ccf <- data.frame(fread(ccf_file))

###########################################################################################

Variant_Type <- c("Missense_Mutation","Nonsense_Mutation","Frame_Shift_Ins","Frame_Shift_Del","In_Frame_Ins","In_Frame_Del","Splice_Site","Nonstop_Mutation")

###########################################################################################
## 关注的病理亚型
if(type == "IGC"){
    dat_info <- subset( dat_info , Type == "IM + IGC" & Type != "IM + IGC + DGC" )
}else if(type == "DGC"){
    dat_info <- subset( dat_info , Type == "IM + DGC" & Type != "IM + IGC + DGC")
}else if(type == "All"){
    dat_info <- subset( dat_info , Type != "IM + IGC + DGC" )
    #dat_info <- dat_info
}

dat_ccf <- subset( dat_ccf , ID %in% dat_info$ID )

###########################################################################################
## 判断基因多少比例出现在Trunk、Pre_Private、Inv_Private
result_plot <- c()
gene_list <- dat_gene$Gene_Symbol

for( geneN in gene_list ){

    tmp <- subset( dat_ccf , Variant_Classification %in% Variant_Type & Hugo_Symbol == geneN )

    if(nrow(tmp) > 0 ){
        tmp$LOH <- ifelse( tmp$minor_cn == 0 & tmp$total_cn == 2 , "LOH" , "Non-LOH" )
        tmp <- subset( tmp , Class != "IM" )

        res_tmp <- c()
        for( sample in unique(tmp$ID) ) {
            tmp_use <- subset( tmp , ID == sample )
            tmp_use <- unique(tmp_use[,c("ID" , "CLS" , "LOH" , "CCF_adj" , "Class")])

            ## 一个人若发生多个突变，算最早的
            if( length(which(tmp_use$CLS=="clonal [share]")) > 0 ){
                tmp_use <- subset( tmp_use , CLS == "clonal [share]" )
            }else if( length(which(tmp_use$CLS=="clonal [early]")) > 0 ){
                tmp_use <- subset( tmp_use , CLS == "clonal [early]" )
            }else if( length(which(tmp_use$CLS=="clonal [NA]")) > 0 ){
                tmp_use <- subset( tmp_use , CLS == "clonal [NA]" )
            }else if( length(which(tmp_use$CLS=="clonal [late]")) > 0 ){
                tmp_use <- subset( tmp_use , CLS == "clonal [late]" )
            }else if( length(which(tmp_use$CLS=="subclonal")) > 0 ){
                tmp_use <- subset( tmp_use , CLS == "subclonal" )
            }
            res_tmp <- rbind(res_tmp , tmp_use)
        }

        res_tmp$CLS_new <- ifelse( res_tmp$CLS %in% c("clonal [share]") , "Share" , "GC branch" )
        res_tmp$Hugo_Symbol <- geneN
        result_plot <- rbind(result_plot , res_tmp)
    }
}

out_name <- paste0(out_path , "/Driver.",type,".CCF.tsv")
write.table( result_plot , out_name , row.names = F , sep = "\t" , quote = F )


###########################################################################################
## 基因分布堆叠图
## 只看存在Share的基因
use_gene <- subset( result_plot , CLS_new=="Share" )$Hugo_Symbol
dat <- subset(result_plot , Hugo_Symbol %in% use_gene )
dat$clonal <- ifelse( dat$CCF_adj > 0.6 , "Clonal" , "Subclonal")

## 计算每个基因克隆和亚克隆的数量
tmp_dat <- dat %>%
group_by( Hugo_Symbol , CLS_new ) %>%
summarize( MutNum = length(which(clonal == "Clonal")))
tmp_dat$CLS <- "Clonal"

tmp_dat2 <- dat %>%
group_by( Hugo_Symbol , CLS_new ) %>%
summarize( MutNum = length(which(clonal == "Subclonal")))
tmp_dat2$CLS <- "Subclonal"

tmp_dat <- rbind( tmp_dat , tmp_dat2 )

## 计算p值
trans <- function(num){
    up <- floor(log10(num))
    down <- round(num*10^(-up),2)
    text <- paste("p == ",down," %*% 10","^",up)
    return(text)
}

result_plot2 <- c()
for( geneN in unique(tmp_dat2$Hugo_Symbol) ){

    tmp_use <- subset( tmp_dat , Hugo_Symbol == geneN )

    a <- subset( tmp_use , CLS == "Clonal" & CLS_new == "Share" )$MutNum
    b <- subset( tmp_use , CLS == "Subclonal" & CLS_new == "Share" )$MutNum
    c <- subset( tmp_use , CLS == "Clonal" & CLS_new == "GC branch" )$MutNum
    d <- subset( tmp_use , CLS == "Subclonal" & CLS_new == "GC branch" )$MutNum

    if(length(a)==0){a=0}
    if(length(b)==0){b=0}
    if(length(c)==0){c=0}
    if(length(d)==0){d=0}

    p <- fisher.test( matrix( c(a,b,c,d) , ncol = 2 ) )$p.value

    if( p < 0.001 ){
        p_text <- trans(p)
    }else{
        p_text <- paste0( "p == " , round(as.numeric(p) , 3) ) 
    }  

    tmp_use <- tmp_use %>%
    group_by( Hugo_Symbol , CLS_new ) %>%
    summarize( MutNum = MutNum , Ratio = MutNum/sum(MutNum) , CLS = CLS )

    tmp_use$p_text <- ""
    tmp_use$p_text[1] <- p_text
    tmp_use$Ratio[is.na(dat$Ratio)] <- 0
    tmp_use$value_text <- paste0( round(tmp_use$Ratio , 2) * 100 , "%")

    result_plot2 <- rbind(result_plot2 , tmp_use )
}

result_plot2$CLS <- paste0( "GC "  , result_plot2$CLS )
result_plot2$CLS_new <- factor( result_plot2$CLS_new , levels = c( "Share" , "GC branch") )
result_plot2$CLS <- factor( result_plot2$CLS , levels = c("GC Clonal" , "GC Subclonal") )

gene_order <- c(
    "TP53" , "ARID1A" , "CDH1" , "APC" , 
    "ERBB2" , "PIK3CA" , "RNF43" , "MAP2K7" ,
    "MTRR" , "MUC6" , "CFTR" , "BMP6" , "GAL3ST3")
result_plot2$Hugo_Symbol <- factor(result_plot2$Hugo_Symbol , levels = gene_order )

plot <- ggplot( data = result_plot2 , aes( x = CLS_new , y = Ratio , fill = CLS ))+
    geom_bar(position = "stack", stat = "identity") + 
    theme_bw()+
    labs(x="",y="Proportion")+
    facet_grid(.~Hugo_Symbol) +
    theme(panel.grid = element_blank())+
    scale_fill_npg() +
    ylim(0,1.05)+
    #geom_text(aes(label=p_text , y = 1.05 ,x = 1.5),parse = TRUE,size=4)+
    geom_text(aes(label=value_text) , position=position_stack(vjust = 0.5) , size=2.5 , color="black")+
    theme(panel.background = element_blank(),#设置背影为白色#清除网格线
                legend.position ='right',
                legend.title = element_blank() ,
                panel.grid.major=element_line(colour=NA),
                plot.title = element_text(size = 12,color="black",face='bold'),
                legend.text = element_text(size = 6,color="black",face='bold'),
                axis.text.y = element_text(size = 7,color="black",face='bold'),
                axis.title.x = element_text(size = 12,color="black",face='bold'),
                axis.title.y = element_text(size = 12,color="black",face='bold'),
                strip.text.x = element_text(size = 7 , face = 'bold'),
                axis.text.x = element_text(size = 8,color="black",face='bold',angle = 45, vjust = 1, hjust=1) ,
                axis.line = element_line(size = 0.5))

out_name <- paste0(out_path , "/Driver_Trunk.",type,".CCF.pdf")
ggsave( out_name , plot , width = 8.5 , height = 3 )


result_plot2 <- subset( result_plot2 , CLS_new == "Share" )
plot <- ggplot( data = result_plot2 , aes( x = Hugo_Symbol , y = Ratio , fill = CLS ))+
    geom_bar(position = "stack", stat = "identity") + 
    theme_bw()+
    labs(x="",y="Proportion")+
    theme(panel.grid = element_blank())+
    scale_fill_npg() +
    ylim(0,1.05)+
    #geom_text(aes(label=p_text , y = 1.05 ,x = 1.5),parse = TRUE,size=4)+
    geom_text(aes(label=value_text) , position=position_stack(vjust = 0.5) , size=2.5 , color="black")+
    theme(panel.background = element_blank(),#设置背影为白色#清除网格线
                legend.position ='right',
                legend.title = element_blank() ,
                panel.grid.major=element_line(colour=NA),
                plot.title = element_text(size = 12,color="black",face='bold'),
                legend.text = element_text(size = 6,color="black",face='bold'),
                axis.text.y = element_text(size = 7,color="black",face='bold'),
                axis.title.x = element_text(size = 12,color="black",face='bold'),
                axis.title.y = element_text(size = 12,color="black",face='bold'),
                strip.text.x = element_text(size = 7 , face = 'bold'),
                axis.text.x = element_text(size = 8,color="black",face='bold',angle = 45, vjust = 1, hjust=1) ,
                axis.line = element_line(size = 0.5))

out_name <- paste0(out_path , "/Driver_Trunk.",type,".Share.CCF.pdf")
ggsave( out_name , plot , width = 6 , height = 3 )
